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Abstract - The collective motion of dust particles during the mode-coupling induced melting of 
a two-dimensional plasma crystal is explored in molecular dynamics simulations. The crystal is 
compressed horizontally by an anisotropic confinement. This compression leads to an asymmetric 
triggering of the mode-coupling instability which is accompanied by alternating chains of in- 
phase and anti-phase oscillating particles. A new order parameter is proposed to quantify the 
synchronization with respect to different directions of the crystal. Depending on the orientation 
of the confinement anisotropy, mode-coupling instability and synchronized motion are observed in 
one or two directions. Notably, the synchronization is found to be direction-dependent. The good 
agreement with experiments suggests that the confinement anisotropy can be used to explain the 
observed synchronization process. 
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Introduction. — Weakly ionized gases containing 
micron-sized dust particles are called complex (dusty) 
plasmas. In the plasma the particles charge up and self- 
arrange enabling formation of strongly coupled and highly 
ordered quasicrystalline phases [1-5] (analogous to colloids 
[6]), called plasma crystals [2,7,8]. In ground-based ex¬ 
periments these crystals are typically composed of plas¬ 
tic microspheres that are injected into a plasma created 
by a radio frequency discharge. The particles charge up 
negatively and levitate in the plasma sheath region above 
the lower electrode where they can form a horizontal two- 
dimensional (2D) monolayer under adequate experimental 
conditions [2,4]. Many dynamical processes can be studied 
rigorously in plasma crystals, in particular, linear [9,10] 
and nonlinear waves [11], resonance effects [12], dynamics 
of dislocations [13-15] and crystal plasticity [16,17]. 

As in many physical, astrophysical and biological sys¬ 
tems [18], cooperative particle motion is an exceptionally 
important element of self-organization in complex plas¬ 
mas. In particular, synchronized motion of particle chains 
was recently discovered in plasma crystals [19]. Synchro¬ 
nization processes in large systems of oscillators have been 
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studied in chemistry, physics and engineering [20], and 
the behavior of chirping crickets [21], or superconduct¬ 
ing Josephson junctions [22] can be described by the Ku- 
ramoto model of globally coupled oscillators [23] which 
can be solved analytically in a mean-field approach. 

In a plasma crystal, the particle-particle interaction is 
strongly influenced by the surrounding plasma. While 
the interaction in the bulk plasma is well described by 
a Yukawa potential [1], the strong ion flow in the plasma 
sheath region distorts the screening cloud [24,25]. This 
plasma wake below the particles adds an attractive com¬ 
ponent to the interaction [26] which was described theo¬ 
retically as a pointlike positive effective charge below each 
particle [27]. Due to the finite vertical confinement of a 2D 
plasma crystal, there is an out-of-plane wave mode which 
has an optical dispersion relation in addition to the two 
in-plane modes with acoustic dispersion. If the vertical 
confinement is smaller than a critical value, the longitu¬ 
dinal in-plane mode and the out-of-plane mode intersect 
and form an unstable hybrid mode in the vicinity of the 
intersection. During this mode-coupling instability (MCI), 
energy is continuously transferred from the flowing ions to 
the crystal, breaking the crystalline order if the damping 
rate is small enough [28,29]. 
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Due to the lattice symmetry, MCI in the shallow cross¬ 
ing regime is equally strong in three directions in a per¬ 
fect hexagonal lattice [29]. In the experiment of ref. [19], 
however, the instability was dominant in one direction. 
Synchronized motion of particle chains was observed. The 
process of synchronization was measured by calculating 
the Shannon entropy of the instantaneous phases of neigh¬ 
boring particles as well as the distribution of frequencies. 
In ref. [19], an inhomogeneity of the horizontal confine¬ 
ment was suggested to be a reason for this asymmetry in 
the crystal, but it was not possible to study the origin of 
the deformation of the crystal in detail. 

The influence of an anisotropy in the horizontal confine¬ 
ment on a rotating plasma crystal was studied in ref. [30]. 
It was shown that even small anisotropies may consider¬ 
ably affect the dynamical behavior of the system. 

In this paper, we demonstrate with simulations that 
an anisotropy of the horizontal confinement can cause an 
asymmetric triggering of MCI. At the onset of the insta¬ 
bility, synchronized particle motion is characterized by a 
new order parameter that is sensitive to the direction of 
the synchronization pattern. Depending on the orienta¬ 
tion of the confinement anisotropy, MCI and synchronized 
motion are observed in one or two directions. 


Experiment. — The experiment of ref. [19] will be 
briefly outlined below. Argon plasma was produced us¬ 
ing a capacitively coupled radio frequency discharge at 
13.56 MHz with a forward power of 12 W. The microparti¬ 
cles formed a monolayer with mean interparticle distance 
a = (480 ± 10) fim. The particle x and y positions were 
obtained with subpixel accuracy from a top-view camera 
operating at 250 frames per second. The axes were chosen 
as depicted in the inset of fig. 1. The gas pressure was 
reduced from 0.94 Pa to 0.92 Pa to initiate the MCI. 

The spectral distribution of particle velocity fluctua¬ 
tions [see eq. (5)] in Aj-space is highly anisotropic [19]. As 
can be seen in fig. 2(c), bright ’hot spots’, the fingerprints 
of the developed MCI, appear in two directions, in contrast 
to a perfect hexagonal crystal where the three directions 
are equally strong [29]. 

Simulation. — Molecular dynamics simulations have 
proven to be an adequate tool to study and compare a wide 
range of experimental conditions [31-34]. The equations 
of motion read 


mvi + muri = ^ Fji ^ Ci^ Li, ( 1 ) 


where Vi is the position of particle i, m the particle mass 
and u the damping rate. 

The force exerted by particle j (and its wake) on parti¬ 
cle i is 


Fji — 




q\Q\ 


( 2 ) 




Fig. 1: (a) Map of the instantaneous phases 0i,6>=o° of the 

particle oscillations at time t = 2.5 s. The particle positions 
are indicated by black dots. The phases of the particles are 
interpolated between the particle positions in order to be visu¬ 
alized in a map. Lines of particles with similar phases appear 
as stripes. The window has a side length 6.5 mm. (b) Map of 
the order parameter [see eq. (6)] at the same time step, 

(c), (d) The same for 0 = 60°. The inset shows the reference 
frame. The direction denoted by angle 0 is shown as a red 
arrow, the line perpendicular to it as a dashed line. 


where Q < 0 is the particle charge, A is the screening 
length, Vji = Vi - Vj and = Vi - {vj - Se^). To 
model the ion wake effect a positive ’extra charge’ g (0 < 
q < IQI) is added a fixed distance S (S < A) below each 
particle. Note that since in general 7^ rwij, the forces 
are nonreciprocal due to the ion wake effect. The ion wake 
(described in detail in [6]) is known to be responsible for 
triggering the MCI [29]. 

To form a monolayer, the equally charged particles have 
to be confined vertically as well as horizontally. In the 
experiment, the confinement can be controlled, e.g., by 
varying the discharge power or gas pressure [29]. In sim¬ 
ulations it is treated as a tunable parameter, allowing 
us to control the crystal stability and anisotropy effects. 
The anisotropic parabolic confinement force in the hor¬ 
izontal plane is characterized by confinement parameter 
f2|| = 27r/|| acting in the direction of angle a (measured 
from the x-axis), and = 27rf_\_ that is perpendicular to 
it. Thus, 

( ft^Xi + ^a{xi COS 2a -h Vi sin 2a)\ 

sin 2a - yi cos 2a) , (3) 

j 

where Vig and Via are the symmetric and asymmetric con¬ 
tributions to the horizontal confinement, and Viz the verti¬ 
cal confinement parameter. The symmetric and asymmet- 
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ric contributions can be expressed as ±r2^)/2. 

The orientation of the confinement anisotropy can thus be 
changed without changing the choice of the axes depicted 
in the inset of fig. 1, leading to a horizontal compression 
of the crystal in the direction denoted by the angle a. 

The particles are also coupled to a Langevin heat bath 
of temperature T = 300 K, 

{Li{t)) =0, {Li{t r)Lj{t)) = 2umTSijS{r). (4) 

Sij is the Kronecker delta and S(t) is the delta function. 

In a simulation run, a system of 16384 particles, each 
with a mass m = 6.1 x 10“^^ kg and charge Q = — 19000e, 
is first equilibrated at /y = f± = 0.145 Hz and a large 
vertical confinement fz = Clzl‘^7r = 23 Hz that prevents 
the onset of MCI. When a crystal is formed in the center 
of the monolayer, the horizontal frequencies are changed 
to /|| = 0.156 Hz and f± = 0.137 Hz to introduce an 
anisotropy. After equilibration, the vertical confinement 
is finally reduced to fz = 20 Hz to trigger the instabil¬ 
ity, this moment corresponds to t = 0. Because of the 
sixfold symmetry of the crystal, it is sufficient to study 
the orientation of the confinement anisotropy in the range 
0° < a < 30°. Here, two simulations with a = 30° and 
(a = 0° are considered. The damping rate is assumed to 
be z/ = 1.26 s“^, the screening length is A = 380//m. A 
pointlike wake charge q = 0.2 jQI is a distance 5 = 0.3 A 
below each particle. 

Analysis methods. — The radial pair correlation 
function in the horizontal plane, ^(r), is used to measure 
the inhomogeneity in the hexagonal lattice. An ellipse is 
fitted to the first six peaks of g{r). The tilt angle [3 and 
the eccentricity e are used to quantify the deformation of 
the crystal. 

The particle current [35] for the longitudinal in-plane 
mode is defined as 

= (5) 

3 

where v^{t) is the component of the velocity of particle j 
at time t parallel to wave vector k = (kx^ky). The par¬ 
ticle current fluctuation spectra of the longitudinal mode 
V{k,f) are then calculated using the Fourier transform. 
To show the spectra in the xy plane, V{k, f) is integrated 
over a frequency range 14 Hz < / < 18 Hz centered on 
the hybrid frequency /hyb = (16 ± 1) Hz. The MCI, where 
the out-of-plane mode couples to the longitudinal in-plane 
mode, appears as hot spots in the spectra of both modes 
[29]. The out-of-plane mode is not considered since it is 
not available for the experimental data. In the simula¬ 
tions, the integrated spectrum of the out-of-plane mode is 
very similar to that of the longitudinal mode. The bor¬ 
der of the first Brillouin zone is calculated from the static 
structure factor S{k) = m where N 

is the number of particles, the sum runs over all pairs of 
particles, and the averaging is performed over time. 



Fig. 2: Synchronized particle motion at the onset of mode¬ 
coupling induced melting of 2D plasma crystal, (a) Pair corre¬ 
lation g(r) in the horizontal plane at t = 0. (b) g(r) in a smaller 
window of side length 1.2 mm with the first peaks shown as 
solid circles. An ellipse (solid line) is fitted to the positions of 
the peaks, its semiaxes of length A = (0.500 ± 0.014) mm and 
B = (0.454 ± 0.011) mm are also shown as solid lines. The tilt 
angle of the ellipse is /3 = (27 ± 2)°. The dashed circle with a 
radius identical to A is shown to guide the eye. (c) Integrated 
particle current fluctuation spectrum of the longitudinal mode 
in the fc/c^-plane, with arbitrary units in a logarithmic scale, 
calculated from the first 3.2 s of the data. The border of the 
first Brillouin zone is shown as a dashed line, (d) Order pa¬ 
rameter Re{t) measuring the degree of synchronization of the 
particles as a function of time t, see eq. (6). 


The chains of synchronized particle motion (see fig. 1) 
cannot be characterized using the Kuramoto order param¬ 
eter = (l/A") [20], because neighboring chains 

tend to be in antiphase. The contributions to the order 
parameter would thus cancel even in the presence of a syn¬ 
chronization pattern. Therefore, we define a local order 
parameter as 


where f>i^o is the phase of the oscillation of particle i in 
the direction denoted by angle 0 at time t and nn is the 
number of nearest neighbors, kj = 0 if particle j is on the 
line passing through particle i perpendicular to direction 
denoted by and kj = 1 otherwise^. In the inset of fig. 1, 

^ Since the crystal is highly ordered, the definition of particle lines 
is straightforward. We consider two neighboring particles i and j to 
be on a line if the angle between the i—j bond and the line is smaller 
than 30°. 


1 f 

RiAt) = — E 

\j=i 
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the direction denoted by ^ = 60° is indicated by an arrow 
and the line perpendicular to it by a dashed line. The 
cosine of the phase differences are thus added for near¬ 
est neighbors on the same line and subtracted for nearest 
neighbors on the subsequent lines, leading to Ri^e = 1 if 
particle i is in a region with perfect alternating in-phase 
and out-of-phase oscillating lines of particles. In the oppo¬ 
site case, Ri^e = —1- If there is no phase relation, Ri^e 0. 

The instantaneous phase (pi^e is calculated from the pro¬ 
jection of position Vi in the horizontal plane onto the di¬ 
rection denoted by 0. The instantaneous deviation from 
the time-averaged particle position is obtained with a slid- 
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ing window of length 0.2 s. The phase is then assumed to 


grow linearly by 27r between each maximum of the devia¬ 


tion. An order parameter for the system is calculated by 


averaging over all particles Re{t) = {Ri,e{t))i. The three 

) 

main directions of the crystal, 0 = 0°, 60° and 120°, are 

V \ A 

considered. 


In ref. [36] a local order parameter was used to increase 
the resolution for a system where the number of oscilla- 




tors is small. The local order parameter proposed here 
is sensitive to the orientation of the synchronization pat¬ 
tern. In fig. 1, maps of <pi^o and Ri^o are shown for the 
experimental data at a characteristic time t = 2.5 s for 
0 = 0° and 0 = 60°. Lines with two different orientations 
are apparent for ^^ 0 = 0 °, see fig. 1(a). The corresponding 
order parameter R^ q^qo [see fig. 1(b)] is sensitive to the 
lines that are oriented along the ^-axis which are located 
in the lower part of the inspection window. For 0 = 60° 
[figs. 1(c) and (d)], the largest values of the order param¬ 
eter are concentrated in the upper part of the window. 


10 15 20 

t[s] 


Fig. 3: Same as fig. 2, but for a molecular dynamics simulation 
of a crystal with an anisotropy in the horizontal parabolic con¬ 
finement. The direction of the largest confinement frequency 
/|l = 0.156 Hz, given by angle a = 30° [see eq. (3)], is indicated 
by arrows in (a) and (c), the frequency in the perpendicular 
direction is f± = 0.137 Hz. The vertical confinement frequency 
is fz = 20 Hz. In (b), the semiaxes of the ellipse fitted to 
the first peaks of g{r) are of length A = (0.498 ± 0.004) mm 
and B = (0.464 ± 0.003) mm. The tilt angle of the ellipse is 
/3 = (29.7 ± 0.5)°. In (c), the first 25 s of the data are used to 
calculate the spectrum. 


Results. — The experimental data of ref. [19] is ana¬ 
lyzed in a window containing about 800 particles near the 
center of the crystal. The pair correlation g{r) is shown in 
fig. 2(a). A deviation from a perfect hexagonal structure 
can clearly be seen. In fig. 2(b), an ellipse is fitted to the 
first peaks of ^(r), its tilt angle is /3 = (27 ± 2)°. The 
value of the eccentricity is e = 0.42 ±0.07. 

The phases <pi^o are calculated for a smaller window of 
side length 6.5 mm containing about 230 particles. In this 
region synchronized particle motion was observed. As can 
be seen in fig. 2(d), the order parameter Rg has significant 
positive values for 0 = 0° and 0 = 60°. In the latter case, 
Rq^qqo increases between t ^ 1 s and t 2 s and then 
saturates at a value of Rq^qqo 0.4. At t ^ 3.2 s the 
crystal melts and the order parameter drops back to zero. 
For 0 = 0°, Re=o° increases much more slowly in a time 
interval 1 s < t < 3 s before also decreasing again when the 
crystal melts. Re=i 20 ° becomes slightly negative during 
the phase of synchronized motion in the other directions 
(see supplementary movie mei-synchronization.mp4 for 
the time evolution of the order parameter). 

In the first simulation, the crystal was compressed at 
an angle oi a = 30°. The region of interest is chosen 
to be of the same size as in the experiments. The pair 
correlation g{r) is shown in fig. 3(a). An ellipse is fitted 


to the first peaks of g{r) [see fig. 3(b)], the value of the 
tilt angle, [3 = (29.7 ± 0.5)°, is close to the experiment, 
the value of the eccentricity e = 0.36 ± 0.03 is slightly 
smaller. The integrated spectrum of the longitudinal mode 
[see fig. 3(c)] shows bright hot spots in two main directions 
of the crystal. The hot spot at = 120° that would be 
expected for a perfect hexagonal lattice is almost absent. 

The main characteristics of the synchronization process 
in the experiment are recovered in the simulations, albeit 
with the roles of ^ = 0° and 0 = 60° interchanged: As 
can be seen in fig. 3(d), the order parameter saturates 
at Re=oo — 0.6 rather quickly in one direction, while it 
follows more slowly in the other. At t 27 s, both order 
parameters decrease. For 0 = 120°, the order parameter 
decreases to about —0.2 to —0.3. Note that the time scale 
is larger than in the experiment. 

The order parameter is slightly negative even at t = 0. 
This could be explained by noting that before the onset of 
MCI, the particle movement is slightly correlated between 
nearest neighbors due to their mutual repulsion. This 
small positive correlation leads to a negative Rq, since 
the phase differences to the four neighbors on the next 
lines are subtracted from the phase differences to only two 
neighbors on the same line. The decrease of Re=i 20 ° dur- 
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ing the period of synchronization can be explained by the 
fact that the three directions of the projection are not or¬ 
thogonal. Consequently, alternating lines of in-phase and 
anti-phase particles in one direction lead to a negative or¬ 
der parameter in the other directions. 

In a second simulation, the orientation of the confine¬ 
ment anisotropy was a = 0°. The pair correlation g{r) 
and the ellipse fitted to the first six peaks of g{r) are 
shown in fig. 4(a) and (b). As expected, the ellipse is 
tilted by only a small angle of (3.0 ± 0.5)°. The eccentric¬ 
ity e = 0.42±0.03 is larger than in the first simulation. The 
integrated particle fluctuation spectrum [fig. 4(c)] shows 
that the MCI is dominant only in the x direction. 

Here, synchronized motion is only observed in ^ = 0° di¬ 
rection, see fig. 4(d). The corresponding order parameter 
increases between t 5 s and t 15 s and subsequently 
saturates at Re=oo — 0.8. Re decreases to negative values 
(about —0.2 to —0.3) in the other two directions 0 = 60° 
and 0 = 120°. 

At t 20 s, Re={)o decreases and increases again at t 
21.5 s. This can be understood as follows. As observed 
in [33], a molecular dynamics simulation of MCI does not 
lead to a complete melting of the crystal but rather to 
cycles of partial melting and recrystallization. Thus the 
synchronization process does not completely stop as in the 
case of the experiment. 

Discussion and Conclusion. — The good agreement 
of the simulation (fig. 3) with the experiment (fig. 2) sug¬ 
gests that the confinement asymmetry can be used to 
explain the observed anisotropic triggering of MCI and 
the synchronization process. The anisotropy of the spec¬ 
tral intensity of the particle velocity fluctuations indi¬ 
cate undoubtedly that the MCI is sensitive to a weak 
anisotropy in the horizontal confinement. The dispersion 
relations for a sheared crystal were examined theoretically 
in refs. [37,38]. 

In the experiment, the hot spot at ^ = 60° is much 
brighter than the one in the opposite direction 0 = 240° 
which almost vanishes, see fig. 2(c). This effect — which is 
much weaker in the simulations — may be due to further 
anisotropies in the crystal structure, stemming for exam¬ 
ple from defect chains near the boundary of the plasma 
crystal. It will be subject to further studies. 

Confining a plasma crystal in the horizontal plane al¬ 
ways makes it internally inhomogeneous. A parabolic 
confinement is often used in the literature, see, e.^., 
[16,17,31,32]. The scaling laws of plasma crystals are also 
well known. In particular, the interaction range k = a/X 
of such clusters is weakly depended on the strength of hor¬ 
izontal confinement parameter 

K oc , (7) 

as is easy to verify by using refs. [15,31,39]. Taking into 
account that the critical vertical confinement for MCI to 
be triggered, ftz,crit^ is known to be strongly dependent on 
the particle interaction range tz (see ref. [29] for details). 


0 200 400 600 800 1000 



X [mm] 



c) 

;:i 

B 

I 

■^-5- 

- 10 - 


m w 

UiX 


W m 


-10 -5 0 5 10 

kx [1/mm] 



Fig. 4: Same as fig. 3, but for a different direction of the 
anisotropy in the horizontal confinement, o = 0°. In (b), the 
semiaxes of the ellipse fitted to the first peaks of g(r) are of 
length A = (0.504 ± 0.004) mm and B = (0.458 ± 0.004) mm. 
The tilt angle of the ellipse is /3 = (3.0 ± 0.5)°. In (c), the first 
20 s of the data are used to calculate the spectrum. 


the influence of the horizontal confinement strength be¬ 
comes apparent. For plasma clusters, the dependence of 
^;z,crit on K is described by [29] 




T(«) 


const, T{k,) = 


± 3k ± 3 


( 8 ) 


The dependence of r^;s,crit on can be calculated by com¬ 
bining eqs. (7) and (8): 


^^2:,crit A / \ 


A{k) 


1 

4 


V k2 + 3k + 3 


(9) 


Since A{k) ^ 1 ai k 1, the relative variation of the insta¬ 
bility threshold is practically proportional to the relative 
variation of the horizontal confinement strength. 

The prediction of eccentricity e due to an anisotropic 
horizontal confinement can be deduced from eq. (7), 
yielding = [l - (K(fl||)/K(fl_L))2] ~ 0.35. Of 

course, in a more detailed analysis the orientation of the 
anisotropy would have to be taken into account, as the 
compressibility of the crystal depends on it. Still, this 
estimate is not far from the values of experiment and sim¬ 
ulations. 

The role of frequency synchronization [19] was not stud¬ 
ied here since special care was taken to quantify the ori¬ 
entations of the phase synchronization processes. The in¬ 
terplay of phase and frequency synchronization during the 
onset of MCI is an important point in the understanding 
of the collective phenomenon. 
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An order parameter motivated by the Kuramoto model 
is often used to quantify synchronization processes [36]. 
If the interaction is repulsive, complex patterns can arise 
that call for a detailed analysis. For example, traveling 
waves [40] or competing domains of different chirality [41] 
were observed. Here, a local order parameter was pro¬ 
posed which is sensitive to the orientation of the observed 
synchronization patterns. 

To conclude, it was shown in simulations that an 
anisotropy of the horizontal confinement can cause an 
asymmetric triggering of MCI which is accompanied by 
particle chains with synchronized motion. To the best 
of our knowledge, it is reported for the first time that a 
horizontal compression in simulations of a plasma crys¬ 
tal reproduces well the synchronization process observed 
in experiments. For an appropriate orientation of the 
anisotropy, MCI is triggered in two directions which leads 
to competing synchronization patterns. If MCI is trig¬ 
gered in one direction, a single pattern dominates. A 
new order parameter was proposed that is able to quan¬ 
tify direction-dependent synchronization. We were thus 
able to identify synchronization patterns that show a pro¬ 
nounced anisotropy. 
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